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Abstract 

We study the production of the four-lepton final state predominantly produced 

by a pair of electroweak Z bosons, ZZ. Using the LoopSim method, we merge NLO QCD 
results for ZZ and ZZ-i-jet and obtain approximate NNLO predictions for ZZ production. 
The exact gluon-fusion loop-squared contribution to the ZZ process is also included. On 
top of that, we add to our merged sample the gluon-fusion ZZ-|-jet contributions from 
the gluon-gluon channel, which is formally of N^LO and provides approximate results at 
NLO for the gluon-fusion mechanism. The predictions are obtained with the VBFNLO 
package and include the leptonic decays of the Z bosons with all off-shell and spin-correlation 
effects, as well as virtual photon contributions. We compare our predictions with existing 
results for the total inclusive cross section at NNLO and find a very good agreement. Then, 
we present results for differential distributions for two experimental setups, one used in 
searches for anomalous triple gauge boson couplings, the other in Higgs analyses in the 
four charged-lepton final state channel. We find that the approximate NNLO corrections 
are large, reaching up to 20% at high transverse momentum of the Z boson or the leading 
lepton, and are not covered by the NLO scale uncertainties. Distributions of the four-lepton 
invariant mass are, however, stable with respect to QCD corrections at this order. 
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1 Introduction 

The production of a pair of electroweak vector bosons constitutes an excellent avenue to test the 
electroweak sector of the Standard Model (SM) at the LHC. This class of processes provides for 
example information on the non-abelian structure of the Lagrangian. Of particular relevance is 
the production of a pair of Z bosons, where the underlying gauge structure of SU{2)l ® H(l)y 
predicts that tri-linear couplings are actually absent at tree-level in the SM. Furthermore, there 
is a contribution from an s-channel Higgs boson resonance, produced in gluon-fusion (GF) via a 
heavy-quark mediated effective Higgs-gluon coupling. This part allows to perform interference 
studies and off-shell Higgs width measurements [lIllElilEllel. 

By the experiments, this process is measured indirectly via Z bosons decaying to a pair of 
charged leptons or neutrinos each. In this article, we will focus on the four charged-lepton final 
state 

PP ^ It(-2 + X. ( 1 ) 

In the past years, the ATLAS and CMS collaborations have provided a rich collection of measure¬ 
ments with increasing accuracy of the four lepton signal, and limits on the trilinear anomalous 
ZZy and ZZZ vertices have been deduced 

The next-to-leading order (NLO) QCD corrections to ZZ production were first computed 
in Refs. mm for on-shell production, and including the leptonic decays and spin correlations 
in Refs. [in US]. They turn out to be sizable, of around 50%, and exhibit a relevant phase- 
space dependence. In differential distributions, much larger corrections up to the order of 10 
can appear. This makes the approximation of correcting the LO differential distribution by the 
total K-factor (the ratio of the NLO over the LO total cross section), an unreliable estimate of 
the true NLO differential cross section and can severely underestimate its size. The origin of 
the large magnitude of the corrections is twofold. At NLO, new partonic sub-processes appear, 
including those with enhanced gluonic parton distribution functions (PDFs). This explains in 
part the size of the total K-factor. On the other hand, some new topologies appear for the 
first time only at NLO. Among them is a topology with a soft or collinear boson emission from 
a quark or anti-quark, which results in an agaEW^T^'^iPTj/'iTT'v) enhancement for a number of 
observables naiiaES]. 

The one-loop gluon-induced corrections gg —>• if ^t^2 + currently known only at 

LO. They were first reported for on-shell production in Refs. m m- Results including the 
leptonic decays naiiH] are also available and studies in the framework of Higgs measurements 
have also been carried out (e.g. Refs. dlElEniEI]). Formally, they contribute to ZZ production 
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only at next-to-next-to-leading order (NNLO) QCD, but due to the large gluon PDFs at the 
LHC, their numerical impact is larger than this naive counting of coupling constants suggests. 
Depending on the selected cuts, their contribution ranges from a few percent up to ten percent. 
This prediction suffers from large scale uncertainties. However, results at NLO QCD for the 
gluonic contributions are expected to be available soon - the real corrections, ZZj production, 
and the virtual two-loop corrections are already known [221 [25] . 

The NLO electroweak corrections for on-shell ZZ production have been computed in Refs. [241 
125] . They yield only a modest contribution, ranging from the few percent level for integrated 
cross sections up to 10-20% for high-pr observables. 

At NNLO QCD, further new partonic channels and new topologies contribute. Thus, to 
match with the expected experimental precision, it is mandatory to assess the size of these 
NNLO corrections, not only at the total cross section level, but also for the differential distri¬ 
butions. ZZj production at NLO QCD provides the mixed real-virtual and the double real 
O(a^) contributions to the NNLO results. They were hrst computed in Ref. [25] for on-shell 
production and account for the new sub-processes and the new topologies appearing for the 
first time at NNLO. Thus, they are expected to provide the dominant contribution to the total 
NNLO prediction for selected observables. 

The NLO QCD corrections to the double real-emission process, pp —)■ (- 2^2 JJ + ^ have 

been reported recently m with corrections around 10%. The size of the corrections is relatively 
mild, if adequate central scales are chosen, due to the absence of new channels and phase-space 
regions opening up at this order, although the uncertainty from varying the factorization and 
renormalization scale gets greatly reduced. 

The two-loop virtual corrections for off-shell ZZ production have been presented in several 
publications [251 ES EHl El] and the NNLO QCD result for the inclusive total cross section has 
been reported recently [52]. The size of the NNLO corrections compared to the NLO result is 
about 10%. Up to date, no fully differential NNLO predictions are available in the literature. 

In this article, we employ the LoopSim method [I51E3] together with the NLO predictions 
for ZZ and ZZj production and the LO ones for GF-ZZ and GF-ZZj, calculated by the Monte 
Carlo program VBFNLO [3411551136] . From this procedure, we obtain merged samples that, 
in certain regions of phase space, are expected to account for the dominant part of the NNLO 
QCD corrections to the ZZ production process. The abovementioned procedure has been used 
recently for other diboson production processes [acs], as well as V-|-jet processes m, resulting 
in corrections ranging from 30-100% for selected distributions. Such sizable corrections should 
be taken into account in experimental analyses. 

The article is organized as follows: In Section]^ the details of the theoretical framework of 
our calculation are given. In section]^ first, we compare our predictions with the existing ones 
presented in Ref. [52] for on-shell Z pair production at the integrated cross section level at NNLO. 
Afterwards, differential distributions are also presented for two set of cuts -- in Section 3.2.1 


the setup of the ATLAS and CMS experimental analyses on ZZ production has been closely 
followed, while in Section |3.2.2[ Higgs search cuts are imposed, following the CMS analysis of 
Ref. |38| . Finally, in Section]^ we present our summary and conclusions. 


2 Theoretical Framework 

Production of the four-lepton final state happens at leading order mainly via the quark-anti¬ 
quark t/u-channel diagram (see Fig. [^(a) for a representative Feynman diagram). The bulk of 
the contribution comes from on-shell ZZ production, as in that case the electroweak coupling 
in the Z decay gets effectively replaced by the corresponding branching ratio, which is about 
10% for a Z boson decaying into a pair of any charged leptons. Instead of Z bosons, also 
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Figure 1: Representative Feynman diagrams for the different production processes contributing 
to pp —>■ li itI 2 ■ Top: qq processes appearing at tree-level. Bottom: GF-initiated processes 
appearing at the one-loop-squared level. 


diagrams appear where these are replaced by virtual photons, denoted as 7 * in the following. 
These typically yield lepton pairs with invariant masses much smaller than the Z mass. Their 
overall contribution strongly depends on the lepton cuts imposed on the final state. Typical 
experimental invariant mass windows for Z bosons have a lower bound of 66 GeV, which reduces 
the 7 * contribution to a negligible level. 

Another possibility to produce this four-lepton final state is production of a single vector 
boson V € (Z,j*), which in turn undergoes a four-body decay into the final-state leptons. An 
example Feynman diagram is depicted again in Fig. ( 6 ). This s-channel contribution is also 
sub-dominant in the SM, since there are no tree-level tri-linear gauge couplings and selection 
cuts on the final-state leptons suppress these contributions due to the limited phase space to 
simultaneously produce the two intermediate vector bosons close to their mass shell. 

Finally, there are one-loop-squared GF diagrams that can also generate the same four-lepton 
final-state. These can either proceed via an s-channel Higgs resonance, which subsequently 
decays into if l^ and is produced via an effective Higgs-gluon-gluon coupling mediated 
by loops of heavy-quarks, predominantly the top quark. The dominant contribution to the 
H —>■ l^li iti-t decay comes from an intermediate ZZ* system, with one on-shell and one 
off-shell intermediate vector boson. The other possibility is a continuum production of two, 
potentially off-shell, Zj^* bosons through a quark-loop box diagram. For typical inclusive cuts, 
the bulk of the GF contributions originates from on-shell ZZ production, similar to the t/u- 
channel diagrams. Production via an s-channel vector boson resonance is forbidden by gauge 
invariance and the Landau-Yang theorem |391 HO] . 

For simplicity, if not stated otherwise, we will refer to the whole process by ZZ production, 
although we will consider all off-shell effects, non-resonant diagrams and spin correlations of 
the four-lepton final state. Also, adding the contributions from virtual photons 7 * is implicitly 
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understood if not mentioned otherwise. 


2.1 ZZ and ZZj production in VBFNLO 

Our calculation relies on the following ingredients. The NLO QCD corrections for the ZZ 
and the ZZj production processes, as well as the LO one-loop gluon-fusion-induced gg —)• 
and gg —>■ contributions are obtained from the VBFNLO package. The 

NLO QCD corrections to ZZj were included first for this study, and, in the following, we give 
some details of the methodology used to compute them to make this work self-contained. We 
follow closely the strategies used for iv'i'j'fj production HU. 

We use the spinor-helicity amplitude method and the effective current approach j32l |33] 
to factorize the electroweak part of the system, containing the leptonic tensor, from the QCD 
amplitude. We first generate the generic qq —)• ViV 2 j amplitudes, where, here and in the 
following, V denotes either a, possibly virtual, Z^*') boson or a virtual photon 7 *. Additionally, 
we also calculate the contribution qq —>■ Vj. Then the leptonic decays Vi —>■ and V —>■ 

are attached via effective currents. Note that, in this way, all the off-shell effects 
and spin correlations are taken into account. All possible flavor and cross-related sub-processes 
are computed from these generic amplitudes. 

At NLO, we need to compute the virtual and the real corrections. To regularize the ul¬ 
traviolet (UV) and infrared (IR) divergences, we use dimensional regularization |44] and the 
anti-commuting prescription of 75 |l5]- We employ the Catani-Seymour algorithm |46j to ex¬ 
plicitly cancel the IR divergences prior to the phase-space integration. 

To evaluate the scalar integrals, we follow the prescription of Refs. HU 1181119] and for the 
one-loop tensor coefficients, we employ the Passarino-Veltman reduction formalism [50| up to 
the box level, but avoiding the explicit appearance of Gram determinants Ha EU, and the 
formalism of Refs. isaisa, with the notation laid out in Ref. (13 pentagon integrals. 

Several checks have been applied to the virtual amplitudes, computed with the package 
described in Ref. Ha) among them, the factorization of the poles, gauge invariance, and 
parametrization invariance Ha- virtual contributions with closed fermion loops, we 

set the mass of all the fermions to zero and check that the poles add to zero. Additionally, we 
have cross-checked both the continuum and Higgs resonance graphs with a second implementa¬ 
tion based on FeynArts and FormCalc [SIESIESIETI- On the amplitude level, there are 10-14 
digits of agreement for each case. The integrated leading order and real emission contributions 
have been checked against Sherpa [58j and agreement at the per mille level has been found. 
Additionally, we have cross-checked that our predictions agree at the per mille level with the 
ones provided in MCFM [SU] . 

A similar strategy is used for the calculation oi ZZ at NLO QCD and the LO gluon fusion 
gg —)■ ifli ^^^2 99 ^2 ^2 9 contributions. We include the leptonic decays via effective 

currents using the spin-helicity formalism and all off-shell effects including Higgs graphs and 
photons are taken into account. More details on the treatment of the challenging numerical 
instabilities appearing in gg —>• ^e found in Ref. [22]. 

For the final-state leptons, we consider Z decays into the first two generations, i.e. electrons 
and muons. We neglect contributions from Pauli interference due to identical particles in the 
final-state, which is below the per mille level. Therefore, all four possible decay combinations 
(Vi, V 2 decaying into e or /r each) give the same result and need to be produced only once at 
the generation level. As we will see later, the cuts are different for same-flavor and different- 
flavor states, and in the Higgs setup case also for electrons and muons, hence after imposing 
analysis cuts, the contributions of the four combinations will differ. For the Z resonance, we 
employ a modihed version of the complex-mass scheme | 6 ^ where the weak mixing angle is 
kept real. We work in the five-flavor scheme and use the MS renormalization of the strong 
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coupling constant, with the top quark decoupled from the running of Ug- Diagrams with a final 
state top quark pair, which would appear in the real emission part of ZZj, are considered a 
separate, experimentally distinguishable process and are therefore discarded. Virtual top-loop 
contributions in contrast are included in our calculation. We consider a massless bottom-quark, 
rrih = 0, except for the closed bottom-quark loop amplitudes, where it is set to its pole mass 
value. The latter is important to correctly account for the negative interference between top- 
and bottom-quark loops in the effective Higgs-gluon coupling, while for the continuum diagrams 
the difference between choosing a massless or massive bottom quark is numerically small. 

2.2 Computing the dominant part of NNLO with LoopSim 

For the calculation of approximate NNLO results for pp —)• if £ 2^2 production, we use the 
LoopSim method [la [ 33 ]. This approach allows us to merge, in a consistent manner, the 
ZZ@NLO and ZZj@NLO samples, provided by VBFNLO, yielding a result which is simulta¬ 
neously NLO accurate both for the ZZ and ZZj production processes. This merged sample 
is expected to provide the dominant part of the NNLO QCD correction to ZZ production in 
phase space regions where additional jet radiation becomes important. 

The NLO ZZj predictions provide the double-real and the mixed real-virtual contributions 
at NNLO of the ZZ result. They are divergent upon integration over the phase space of the 
real partons and need to be supplemented with double-virtual contributions to yield the finite 
NNLO result for ZZ production. LoopSim constructs approximate versions of such double- 
virtual terms by utilizing the fact that their structure of divergences has to match exactly that 
of the higher multiplicity contributions. This procedure guarantees finiteness of the combined 
result, while providing a dominant part of the NNLO result for a number of relevant observables. 
In the process of determining the exact divergent terms of the two-loop corrections, some finite 
pieces are generated, which are not guaranteed to match the exact constant part of the two- 
loop diagram. They are, however, proportional to the LO born kinematics and are therefore 
negligible for distributions which receive sizable corrections at NLO. 

To link the two programs, we have made use of an interface [Ej which consists, on one hand, 
of an extension of the VBFNLO program that writes down events in the Les Houches event 
(LHE) [E] format at NLO and, on the other hand, of a class in the LoopSim library that reads 
and processes the LHE events. 

The LoopSim method proceeds in the following steps. First, an underlying structure for 
each NLO ZZj event is determined with the help of the Cambridge/Aachen (C/A) |62l [63] 
algorithm, as implemented in Fast Jet [GllISS], with a certain radius iiLS- This step establishes 
the sequence of emissions in the input event. For simplicity, we combine each pair of oppositely- 
charged leptons to a virtual vector boson V and LoopSim processes diagrams at this level. 
Thereby, we use information from the event generation step to identify the leptons connected 
by a continuous fermion line. For s-channel-type events, which are dominated by contributions 
where only a single boson is attached to the QCD part of the amplitude, one might consider 
to also reflect this in the LoopSim part and combine the four leptons into a single particle. 
The question would then be how to define the transition between the two regions. As this 
contribution is strongly suppressed by the cuts applied later, we do not pursue this further, but 
instead always combine the four leptons into two Z bosons. In the next step, the underlying 
hard structure of the event is determined by working through the ij —>■ k recombinations in 
order of decreasing hardness, defined by the kt algorithm measure |66l 167] . The first nj, particles 
associated with the hardest merging are marked as “Born”. The number of Born particles is 
fixed by the number of outgoing particles in the LO event and is equal to two for the case 
of ZZ production. At NNLO, the Born particles can be either both vector bosons, a boson 
and a parton, or two partons. The remaining particles, which are not marked as “Born”, are 
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then “looped” by finding all possible ways of recombining them with the emitters. This step 
generates approximate one- and two-loop diagrams. 


In the next step, a double counting between the approximate one-loop events generated by 
LoopSim and the exact one-loop events coming from the NLO sample with lower multiplic¬ 
ity is removed. This is done by generating the one-loop diagrams from the tree level events 
first, and then using them to generate all possible one and two-loop events. This set is then 
subtracted from the full result, which amounts to removing the two-loop diagrams with both 
loops simulated by LoopSim and leaving only those with one exact and one simulated loop. To 
satisfy unitarity, these diagrams are assigned a weight equal to the original event times a pref¬ 
actor (_i)™mber of loops^ This guarantees that the sum of the weights of the LoopSim-generated 
events is zero m- For a fully inclusive observable, the integrated cross section of our approxi¬ 
mated NNLO result would be equal to the NLO one (hence the pure contributions would 

vanish). However, in realistic situations, with finite fiducial volumes or in the case of differen¬ 
tial distributions, some of the events generated by LoopSim are removed by cuts or reshuffled 
within histograms, which results in non-vanishing correction and that leads to a genuine, 

approximate NNLO correction. 


The jet C/A and kt algorithms mentioned above depend on the radius Rls, which is a 
parameter of the LoopSim method. The smaller the value of RlS: the more likely the particles 
are recombined with the beam, the larger RhS^ the more likely they are recombined together. 
The value of i?LS is irrelevant for collinear and soft radiation. It affects only the wide angle (or 
hard) emissions where the mergings between particles i and j compete with mergings with the 
beam. In our study, we shall use i?LS = 1) and we shall vary it by ±0.5. The i?LS uncertainty 
will therefore account for the part of the LoopSim method which is related to attributing the 
emission sequence and the underlying hard structure of the events. 


In order to distinguish our predictions with simulated loops from those with exact loop 
diagrams, we denote the approximate loops by n, as opposed to N used for the exact ones. 
With that notation, for processes whose contributions start at tree level like qq ^ ZZ, hLO 
denotes the correction with simulated one-loop diagrams, and hNLO is a result with exact 
one-loop and simulated two-loop contributions. However, for processes that start contributing 
only at one-loop, like gg —)• ZZ, hLO denotes the correction with respect to that first, non¬ 
trivial result. Hence, it is formally an N^LO contribution with respect to the full process of ZZ 
production. 


The GF contribution formally first contributes at NNLO, and consequently we also include 
it in our merged hNLO sample generated with LoopSim. Due to the large gluonic PDFs, this 
process can contribute relevantly despite the suppression. Hence, and since it is gauge 
invariant on its own, by now a common approach in the literature is to add this contribution 
already to the NLO results. We follow this convention, but make the addition explicit by using 
the label “NLO±LO-GF” in this case. Additionally, as mentioned before, we also merge the 
real radiation process GF-ZZj computed at LO to the GF-ZZ result, yielding a contribution 
appearing only at N^LO. Our results will in general also include this contribution, where we label 
the full results as “hNLO ±hLO-GF”. In summary, the GF contribution is always implicitely 
understood to be included at the corresponding order given by the power counting of coupling 
constants, in particular LO-GF in the hNLO result. If additional GF contributions are added, 
this is made explicit in the label. 
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o-LO [pb] 

5.0673(4) +1'^ {Ref. JEM- 5-060 

o-NLO [pb] 

7.3788(10) {Ref. 7.369 

C^NLO+LO-GF [pb] 

7.946(3) 

O'NNLO [pb] 

{Ref. : 8.284 +|°^“) 

0-nNLO [pb] 

8.103(5) (f) -iS (-Rw) 

C^nNUO-l-nLO-GF [pb] 

8.118(5) +S i^^) {Rls) 


Table 1: Comparison with Ref. |32] of total cross sections for on-shell ZZ production at the 
LHC running at \fs = 8 TeV. The errors in brackets are the statistical error from Monte Carlo 
integration, while the percentages give the scale variation error, obtained from changing [ip and 
\xp independently within the range where the ratio hf/^ir is constrained to stay 

within [|;2]. For the hNLO results we additionally give the error due to a variation of Rls 
between 0.5 and 1.5. 

3 Numerical Results 

3.1 Comparison with inclusive NNLO calculation 

We start by comparing the results obtained with LoopSim+VBFNLO with the calculation of 
the inclusive NNLO cross section of Ref. |32j . Here, in contrast to the rest of the paper, we 
use the settings as those of Ref. [32j, i.e. the Z bosons are on-shell and do not decay, hence 
no cut is placed on the final state. Also, all numerical values of masses and couplings are 
taken from there, namely mw = 80.399 GeV, mz = 91.1876 GeV, Gp = 1.16639 x 10“® GeV“^, 
mt = 173.2 GeV, mn = 125 GeV and fJ-R = Up = mz- As PDFs, the MSTW 2008 set [ 68 ] is 
chosen, evaluated at each corresponding order. Since our matrix elements include all off-shell 
effects and spin correlations in our prediction of production, for this comparison, we 

have to modify our code. All the s-channel contributions, like the example diagram in Fig. [^(5), 
or the t/u-channel contributions with 7 * in the intermediate state are set to zero. Hence, the 
only contribution comes from t/u-channel diagrams of Fig. [^with on-shell Z bosons. For the 
GF part, both Higgs and continuum diagrams contribute, but for the latter we also have to 
remove all diagrams with virtual photons. In the phase-space generator, the (Breit-Wigner) 
distributions for the invariant mass of the Z bosons are replaced by ^-distributions at the Z 
pole mass. As leptonic decays of the Z bosons are still simulated internally, we finally need 
to account for this by dividing the result by BR(Z —>■ . The resulting cross sections at 

8 TeV are shown in Table [TJ 

At LO and NLO, agreement at the per mille level is found. The NLO integrated K-factor, 
defined as the ratio of NLO/LO predictions, is 1.46. For comparison, we also show the result of 
adding the GF contribution, formally NNLO, to the NLO result, evaluating both with NNLO 
PDFs. This is the best currently available estimate without requiring the evaluation of two- 
loop diagrams or merging different jet multiplicities. We see that these give an additional 7.7% 
contribution to the NLO cross section, or a total K factor of 1.57 comparing NLO-I-GF to the LO 
result. Additionally, they increase the scale variation uncertainty significantly. The latter comes 
from the fact that the GF contribution is the lowest order accuracy result for the gg channel. 

With respect to NLO, the total NNLO correction, computed in Ref. [H^ and quoted in the 
fourth row of Table is about 12%. Hence, the GF contribution provides the leading part of 
those, namely 60%. Our approximate hNLO and hNLO -|-hLO-GF results are shown in the 
last two rows of Table Their overall agreement with the full NNLO result is good with the 
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difference at the level of 2% only. This is not a priori guaranteed by our method and is consistent 
with the assumption of a LO x effect coming from the genuine finite pieces of the exact two- 
loop virtual amplitudes. These terms are not properly determined by our method, but they are 
covered by the remaining scale uncertainty. The scale uncertainty of our nNLO result is similar 
to NLO-I-GF, and not reduced like for the NNLO result, since the LoopSim method does not 
attempt to reconstruct higher-order terms proportional to the scale dependence in order not to 
underestimate the variation, although, technically, this would be possible. 


3.2 Differential distributions 


In the following sections, results at the LHC at ^/s = 8 TeV will be given for two different sets 
of cuts. In Section 3.2.1, we closely follow the ATLAS and CMS experimental analyses on ZZ 
production, while in Section 3.2.2, we impose Higgs search cuts, following the CMS analysis of 
Ref. |38j . Below, we describe the common settings. 

As input parameters, we use 


mz = 91.1876 CeV , Gp = 1.16637 x 10"® CeV^ , 

mw = 80.398 CeV, a~^ = 132.3407, 

ruH = 125 CeV , sin2(6liy) = 0.22265 , (2) 

Gz = 2.508 CeV , Fp = 0.004017 CeV . 


The mass of the top and bottom quarks, which run in the closed fermion loops, are set to 

mt = 172.4 CeV , nib = 4.855 CeV . (3) 

All other quarks, including external bottom quarks, are taken as massless. 

The jets are defined with the anti-fct algorithm |69| . as implemented in Fast Jet IMl EO], 
with the radius R = 0.4. Independently of the order of a prediction, we use the NNLO 
MSTW2008 [68] PDF set, provided by the LHAPDF [7T| implementation with as{mz) = 
0.11707. 

At fixed order in perturbation theory, the cross section depends on the renormalization and 
factorization scale. As central values for both of those scales, we choose the scalar sum of the 
transverse energy of the system 

f^F,R — ^^0 — 2 ^ Pr.partons + '\JPt,Vi ^Vi 4” Pt,V 2 4” ^^ 2 ) ’ 

where Pt,Vi 2 ^^*4 niy^ 2 ^^le transverse momenta and invariant masses of the recombined, 
opposite-signed charged lepton pairs, respectively. The scale uncertainty is obtained by varying 
simultaneously the factorization and renormalization scale by a factor two around the central 
scale. Additionally, to assess the uncertainties associated with the recombination method used 
by LoopSim, we show the uncertainty bands associated with variations of ±0.5 around Rls = 1- 


3.2.1 ZZ analysis 


In the analysis ol ZZ production, we use the following cuts, inspired largely by the ATLAS 
paper [7|. The settings of the corresponding CMS analysis |8| are comparable. The transverse 
momenta and pseudorapidities of leptons, as well as those of jets (for observables exclusive with 
respect to jet activity), and the distance between leptons and leptons and jets are required to 
stay in the following fiducial volume: 


pt^i > 20 CeV , 
ptjet > 25CeV, 
> 0.3, 
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|%| <2.5, 

Ihjetl < 4.5 , 
AR^ £ > 0.2 . 


(5) 







ZZ ZZ* 

(^LO [fb] 

9.394(9) +2-2% 1.0134(16) 

<7nlo [fb] 

12.057(19) 1.314(3) tlH 

C^NLO+LO-GF [fb] 

12.929(19) tin 1-365(3) tlS 

C^nNLO [fb] 

13.15(8) tin (-“) -aei 1.417(12) + 2 - 0 % + 0 . 8 % 

<7nNLO-l-fiLO-GF [fb] 

13.15(8) (^) +0.9% 1.427(12) +2-3% +o.9% 


Table 2: Inclusive cross sections at ^/s = 8 TeV for the process pp —)■ £^£^£^£2 using the 
cuts of Eq. ([^, separated into ZZ and ZZ* event categories defined in Eq. Q. The errors in 
brackets are the statistical error from Monte Carlo integration, while the percentages give the 
scale variation error, obtained from varying p = PF = ^ [ 5 /^ 0 ; “^Po], with po given by Eq. 

Eor the hNLO results, we additionally give the error due to a variation of Rls between 0.5 and 
1.5. 


To reconstruct the Z bosons from the leptons, we employ the following algorithm. Eirst, all 
invariant-mass pairs of same-flavor and opposite-sign lepton pairs are formed. If all leptons are 
of the same generation, there are in total four possibilities, while for different generations only 
two exist. The invariant-mass pair closest to the physical Z boson mass is labeled Zi and it is 
required to satisfy the cut 66 GeV < minv,Zi < 116 GeV, otherwise the event is discarded. If 
the second pair of leptons, which we denote as Z 2 , falls into the same mass window, the event 
is labelled as ZZ, if not, it is called a ZZ* event, provided that mz 2 > 20 GeV. If the latter is 
not satisfied, the event is rejected. Hence, our two selection types can be summarized as 


ZZ selection: mzi , rnz2 ^ ( 66 , 116 ) GeV , 

ZZ* selection: mzi G ( 66 , 116 ) GeV , mz 2 £ (20, 66 ) U ( 166 , mz, max) GeV , 


( 6 ) 


where mz,max is the maximal mass that can be obtained for a given energy of the system of the 
incoming partons. Note that the terminology adopted for our study differs slightly from that of 
Ref. [7], where ZZ* was used for the union of both categories defined in Eq. Q. 

In Table we present the inclusive cross section at different levels of accuracy. As one 
can see, the overall behaviour is similar to what we have already observed for total on-shell 
ZZ production in the previous subsection. The GE contribution gives a significant correction 
to the NLO result of about -|-7.2% for the ZZ case, while for ZZ* it is only -|-3.9%. In both 
cases the scale dependence is strongly increased. The additional integrated hNLO corrections 
are modest with 1.7% and 3.8% for the ZZ and ZZ* cases, respectively. The dependence on 
the LoopSim-Parameter iiLs is clearly smaller than the remaining scale variation error. 

A more important aspect for our method are, however, differential distributions, where the 
effects can be much larger. To this, we will turn next. 

Fig. [2] shows distributions of the effective mass observable, Ht, defined as a scalar sum of 
the transverse momenta of leptons and jets, and missing transverse energy 

Hf = E PT,jets + PT,1 + Et, miss ■ (7) 

The left panel corresponds to the ZZ and the right panel to the ZZ* types of cuts. In the 
former case, the K factor is very large, both at NLO and at hNLO. In the case of ZZ* selection, 
the K factor is visibly smaller. The leading correction to this observable at NLO comes from 
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Figure 2: Differential cross sections and K factors for the effective mass observable Ht, defined 
in Eq. Q, for the LHC at y/s = 8TeV. The left and the right plot correspond to the ZZ and 
ZZ* selections defined in Eq. (©. The bands correspond to varying fip = by factors 1/2 
and 2 around the central value from Eq. Q. The cyan solid bands give the uncertainty related 
to the i?LS parameter varied between 0.5 and 1.5. The distribution is a sum of contributions 
from same-flavor decay channels (4e and 4/r) and the different-flavor channel (2e2/i). 



Figure 3: Example diagrams contributing to ZZ production at LO, NLO and NNLO. 


configurations, shown in the middle diagram of Fig. with one of the bosons emitted collinearly 
and the other with large transverse momentum recoiling against a hard jet Prjet — Pt,z m, 
which results in a dependence given by 


da- PT,jet 

— oc In-^ 

aii 


( 8 ) 


A similar enhancement occurs at NNLO, where both Z bosons are allowed to be soft or collinear 
and the result is dominated by the dijet type configurations shown in Fig.(right). This explains 
both why the K factors grow with transverse momentum and shows that the rate of this growth 
depends on the selection of the vector boson mass. 
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Figure 4: Differential cross sections and K factors for the effective mass observable Ht, defined 
in Eq. Q for the LHC at ^/s = 8TeV. The left plot corresponds to the part of ZZ* selection, 
Eq. Q, where 20 < mz < 66 GeV, and the right plot to the part of ZZ* selection with 
mz > 116 GeV. Other details are as in Fig. 


In Fig. we split the result of Fig.J^ (right) into the two separate mass regions for the off- 
shell bosons: 20 < mz^ < 66 GeV Fig. |^(left) and > 116 GeV Fig. (right). We see that 
the NLO and the hNLO K factors are much larger in the former case, as putting a smaller mass 
in Eq. Q leads to a stronger logarithmic enhancement. We also see that, at large Ht, it is the 
right plot of Fig. [^that contributes more to the sum shown in Fig. (right), in terms of absolute 
values. This comes from the fact that the born cross section is already significantly larger in 
this case. Therefore, Fig. (right), with larger mz values, which based on Eq. Q results in 
a lower K factor, dominates the total result of Eig. (right). That is why the enhancement is 
smaller there, as compared to Fig. (left), where only the mass region around the Z peak is 
included. 

Let us now turn to Fig. where we show distributions of the transverse momentum of the 
leading Z boson and that of the leading lepton for the ZZ selection. We see that, in both cases, 
the nNLO correction is significant, reaching up to 20% with respect to the NLO result at high pT- 
This magnitude of nNLO correction at high pT is similar to the one already found in WZ |14j 
and WW [IS] production. We note that, above 200 GeV, the R^s uncertainty is much smaller 
than the uncertainty coming from variations of the factorization and renormalization scales. 
Hence, the error related to assigning emission sequence by LoopSim is negligible. The overall 
theoretical uncertainty decreases by 20-30% as we go from NLO to nNLO for these distributions. 
We consider the nNLO predictions shown in Fig. as one of the highlights of our study. We 
emphasize that the inclusion of QGD corrections of that size should be mandatory in all related 
diboson analyses at high transverse momenta. Such large corrections are of huge importance and 
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Figure 5: Differential cross sections and K factors for the transverse momentum of the leading 
Z boson (left) and the leading lepton (right) for the ZZ type selection of Eq. (© at the LHC 
with \fs = 8TeV. All details as in Fig. 


should be accounted for both in computations of SM backgrounds and in searches for anomalous 
couplings. 

We notice that the kink around 900 GeV in the distribution is due to the cut of 

Eq. (§. For a high-pT boson, the decay products are highly collimated with the Z transverse 
momentum and they share approximately equal amounts of its pT- Hence, the mass of the 
dilepton system can be approximated by ~ \p\ ^ARj The imposed phase-space cut 
ARi^i > 0.2 then leads to the condition pT,z ^ 10Since the mi^£ distribution is strongly 
peaked at the Z boson mass, we obtain that the typical separation between the leptons drops 
below 0.2 at pT,z — 900 GeV. This leads to the kink observed in the LO distribution of Fig. 
Also, because the two Z bosons are back-to-back in LO conhgurations and hence have the 
same pt, the same effect happens simultaneously for both Z bosons. With additional parton 
radiation, this effect is smoothed out. The sub-leading Z boson will in general have a smaller 
Pt, as some transverse momentum is carried by the parton. The kink at about 800 GeV in the 
distribution, shown in the right panel of Fig. [51 is also due to the ARi^n cut. We have 
checked explicitly that changing the ARg^g cut to O.T moves the kink to 450 and 400 GeV for 
the and distribution, respectively, consistent with the above discussion. 

Another interesting distribution is the mass of the ZZ system or, equivalently, the mass 
of the system of four leptons, m 4 i, which is shown in Fig. This observable is particularly 
important from the point of view of anomalous triple gauge coupling (aTGG) searches, as the 
effects of anomalous couplings are enhanced in events where large momentum is transferred 
through the triple-boson vertex mu- Such events result in a large invariant mass of the di¬ 
boson system. It is interesting to find that the m^i distribution receives only modest nNLO 
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Figure 6: Differential cross sections and K factors for invariant mass (left) and the transverse 
momentum (right) of the four-lepton system for the ZZ -type selection of Eq. Q at the LHC 
with ^/s = 8TeV. All details as in Fig. 


corrections from QCD, which stay always below 5%. This is to be compared with typical sizes 
of electro-weak corrections, which are not accounted for, and can be of the order of 10% in the 
tail of the distributions, and with the PDF uncertainties, estimated at 5-10%. Hence, compared 
to the above sources of uncertainties, the nNLO QCD corrections are small and we conclude 
that the distribution becomes stable at this order and can be safely used for setting aTGC 
limits [71|8]. 

Fig.[§ (right) shows the distribution of the transverse momentum of the four-lepton system. 
At LO, this distribution is just a (5-function at pr, 4 £ = 0, since there is no jet the four-lepton 
system could recoil against. We do not show this first bin in the figure. The first non-trivial order 
for this observable is then NLO. As we go to nNLO, the distribution receives significant 
corrections of the order of 60-80% for the range presented in the plot, consistent with the NLO 
predictions shown for ZZ-|-jet production in Ref. [26]. 

Finally, in Fig. we show the distribution of the xjet variable defined as 


Xjet 


Sfcsfjets} ^T,k 

^^fcSfjetSjZs} ^T,k 


(9) 


which has been introduced in Ref. m in the context of a dynamical jet veto. Large values of 
this variable correspond to configurations where most of the energy of the final state is carried 
by jets. For the case of ZZ production, a:jet = 0 at LO and the distribution starts to be non¬ 
trivial at NLO. Here we now also show the hrst bin, which almost exclusively consists of the 
xjet = 0 contributions, and which we have scaled down by a factor 30 to fit the plot range. At 
LO, including LO-GF, the result is simply the inclusive cross section of the process. At higher 
orders, this bin corresponds effectively to the contribution after placing a jet veto on the process. 

As we see in Fig.[^ at NLO, where one jet is allowed, the xjet distribution is peaked around 
0.1-0.15, hence, for most events, the jet carries 10-15% fraction of the final state energy. We 
also see, however, that there is a non-negligible tail reaching out to xjet = 0.5. When we move 
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Figure 7: Differential cross sections and K factors at various perturbative orders for the xjet 
variable defined in Eq. (§. The results correspond to the LHC with ^/s = 8TeV and ZZ -type 
selection of Eq. All details as in Eig. 


to nNLO, the yield of xjet events increases over the entire range of the distribution at the 
expense of the energy carried by the Z bosons. In particular, the tail receives corrections with 
K factors of the order of 5 around xjet = 0.5. The nNLO result suggests that the xjet cut, 
used in dynamical jet veto analyses [73], should be set around xjet = 0.2, somewhat lower than 
what could be inferred from the NLO distribution. The i?LS uncertainty is negligible for this 
observable and the renormalization and factorization scale uncertainty is comparable at NLO 
and nNLO. 

3.2.2 Higgs analysis 

We turn now to the discussion of the four-lepton production in the context of Higgs analyses. As 
explained in Section]^ and shown in Fig. the same ZZ —)• dt' signature can be produced with 
and without the intermediate Higgs boson. The latter constitutes an irreducible background 
from the point of view of studies of Higgs production in gluon-gluon fusion, therefore, its precise 
determination is of utmost importance. The analyses optimized for Higgs studies use slightly 
different set of cuts with respect to those employed for aTGC searches. 

Following the CMS analyses of Ref. m, where Higgs properties are studied, we require 

pt,e > 7GeV, \r]e\ < 2.5, 

Pf,^>5GeV, |r/^|<2.4, (10) 

Pt.4arde.t > 20 GeV , mu > 100 GeV , 

7*ii4econd-hardest ^ ^0 GeV , 

and 

40 < mu <120 GeV for the pair with mass closer to mz, 

12 < mu < 120 GeV for the other U pair, (11) 

mu > 4: GeV for any oppositely-charged pair of leptons. 
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Figure 8: Differential cross sections and K factors for transverse momentum of the leading 
lepton (left) and the leading Z boson (right) for the Higgs-type selection defined in Eqs. ( [lOl ) 
and (11). The results correspond to the LHC at ^/s = 8TeV. All other details as in Fig. 


We start by showing in Fig. the distributions of transverse momenta of the leading Z 
and the leading lepton, focusing only on the GF part. The LO result corresponds exactly to 
the diagrams of Fig. (c) and (d). The hLO correction is computed with LoopSim using the 
gg 4£ + j result of Ref. [22] provided by VBFNLO. We see that the correction to and 
is large practically over the entire range shown in Fig. The K factor reaches up to 50% 
for central values and the scale variation does not decrease as we go from LO to hLO (bottom 
panel). It is important to notice that the Rls uncertainty is negligible for those distributions, 
hence our prediction is very reliable within the uncertainty coming from the factorization and 
renormalization scale. 

Let us now turn to the distribution of the invariant mass of the ZZ pair. As argued in 
Ref. [T], even though the Higgs width is extremely small in the Standard Model, interference 
effects between continuum and Higgs-mediated contributions (diagrams (c) and (d) in Fig. 
lead to enhanced four-lepton mass spectra in the region m^i > 2mz- Hence, the off-shell effects 
cannot be neglected, and, as shown recently ms], by comparing the yield at the Higgs mass 
peak with that off the peak, they can be used for setting bounds on the Higgs decay width. 

Our predictions for the four-lepton mass spectra are shown in Fig. Similarly to the ZZ 
selection discussed in the previous section, also here, nNLO QCD effects bring a very small 
correction to this observable. Hence, the 17141 distribution is stable at this order. This is 
important from the point of view of the off-shell effect studies in GF fusion, since a precise 
determination of the expected theoretical yield has impact on setting the Higgs width limits. 
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9: Differential cross sections and K factors for the invariant mass of the four-lepton 
for the Higgs-type selection defined in Eqs. (10) and ©• The results correspond to the 
^/s = 8TeV. All other details as in Fig. 


4 Summary 

We have studied the ZZ production process at the LHC beyond NLO in QCD by merging the 
ZZ and ZZ+jet NLO samples with help of the LoopSim method. We have included the exact, 
loop-induced GF predictions, which are part of NNLO, as well as the GF ZZ-|-jet contribution, 
which is formally of N^LO. The NLO samples were obtained with the VBFNLO package. The 
leptonic decays of the Z bosons, including all off-shell and spin correlation effects, have been 
fully taken into account. 

We have compared the fully inclusive cross section predictions from our framework with the 
exact results computed in Ref. [32j and found a very good agreement within 2%, covered by the 
remaining scale uncertainties, despite the fact that the LoopSim method is missing some finite 
parts originating from the two-loop virtual contributions. Following closely two experimental 
setups, one used in the SM ZZ production and aTGC searches, and the other used in Higgs 
analyses, we obtained results for a selection of differential distributions. 

For the observables sensitive to QCD radiation, the corrections exceed the errors bars of 
the NLO predictions combined with the LO-GF results and range from 20%, in the case of the 
transverse momentum of the Z boson or the leading lepton, up to 100% for the effective mass 
variable Ht- A study of the xjet observable, defined as ratio of the transverse energy sum of the 
jets over the sum of transverse energies of jets and Z bosons, shows that the nNLO corrections 
start to increase quickly when this variable exceeds 0.2. Therefore, when this variable is used 
to impose a dynamical jet veto, a cut should be placed at around this value. For Higgs-type 
selection, we investigated also the ZZ production via GF at nLO. The size of the corrections 
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are around 50% for the transverse momenta of the leading Z and 20% for those of the leading 
lepton. 

For observables which favor the LO kinematics, like m^i, the approximated NNLO QCD 
corrections are small, of the order of 5%, and comparable with the size of the remaining scale 
or PDF uncertainties. 

Modifications to the VBFNLO program used in this article are available on request and will 
be part of a future release. The LoopSim library, together with the Les Houches Event interface, 
is publicly available at https://loopsim.hepforge.org, 
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